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Abstract 


In  many  complex  flows,  large-eddy  simulation  is  difficult  due  to  the  simulta¬ 
neous  presence  of  a  variety  of  flow  features,  often  with  quite  different  resolution 
requirements.  For  example,  a  typical  flow  over  an  airfoil  near  maximum  lift  in¬ 
cludes  laminar,  transitional,  and  turbulent  boundary  layers,  flow  separation, 
unstable  free  shear  layers  and  a  wake.  In  such  situations,  unstructured-grid 
methods  can  gain  great  efficiencies  over  structured-grid  methods  (factor  of  27) 
by  placing  points  based  on  local  resolution  requirements,  rather  than  along  lines 
(as  is  required  with  structured-grid  methods).  For  this  reason,  unstructured- 
grid  large-eddy  simulation  techniques  are  being  developed  using  the  finite  ele¬ 
ment  method  to  solve  the  compressible  Navier-Stokes  equations  with  a  dynamic 
model  of  the  subgrid-scale  stresses.  As  the  purpose  of  this  report  was  to  sup¬ 
port  graduate  student  education,  the  bulk  of  this  report  is  taken  from  the  work 
of  the  student  supported  by  this  research  (Andres  Tejada-Martinez)  .  That 
student  has  focused  on  the  characterization  of  filters  used  within  the  dynamic 
subgrid-scale  model  for  large  eddy  simulation.  Though  he  is  not  quite  finished, 
this  report  clearly  shows  that  much  was  learned  under  the  three  years  of  support 
by  the  AFOSR.  This  grant  also  allowed  partial  support  for  Michael  Yaworski. 
He  is  currently  finishing  his  masters  in  RANS  simulations  and  will  go  on  to 
link  Andres  work  to  combined  LES/RANS  models.  We  have  decided  to  present 
only  Andres  work  here  in  the  interest  of  completeness  and  compactness. 
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1  Introduction 


Large-eddy  simulation  (LES)  is  a  technique  for  computation  of  turbulent  flows 
where  the  large-scale  component  of  the  flow,  carrying  most  of  the  energy,  is 
resolved  by  the  computational  method,  and  the  small-scale  field  motions  are 
modeled.  Previous  to  LES,  the  two  major  approaches  to  simulating  turbulent 
flows  were  Direct  Numerical  Simulation  (DNS)  and  Reynolds- Averaged  Navier- 
Stokes  Simulation  (RANSS).  In  DNS,  the  computational  method  resolves  all  of 
the  turbulent  motions  from  the  largest  scale  down  to  the  scale  where  motion 
is  converted  to  heat  via  viscous  dissipation.  It  is  well  known  that  DNS  is  lim¬ 
ited  to  low  Reynolds  number  flows  due  to  the  increasing  range  of  small  scales 
with  increasing  Reynolds  number.  In  RANSS,  the  space-time  computational 
grid  is  too  coarse  to  resolve  the  flow  instabilities  that  lead  to  and  characterize 
turbulence.  The  grid  is  only  fine  enough  to  resolve  the  mean  flow,  thereby 
requiring  some  type  of  modeling  of  the  statistical  effects  of  all  turbulent  fluctu¬ 
ations  on  the  mean  flow.  LES  is  a  compromise  between  DNS  and  RANSS.  The 
computational  grid  is  sufficiently  fine  to  resolve  some  flow  instabilities,  but  not 
fine  enough  to  resolve  the  energy-dissipating  motions.  The  idea  behind  LES 
is  that  motions  that  are  resolved  are  the  important  ones  and  the  errors  in¬ 
duced  by  modeling  the  small-scale  motions  are  significantly  smaller  than  those 
incurred  in  RANSS.  The  constant  coefficient  Smagorinsky  model  (see  [15]), 
analogous  to  the  linear  viscosity  model  in  RANSS,  is  often  used  to  account  for 
the  un-resolved  small-scale  motions.  A  more  accurate  procedure  is  often  used 
in  which  the  constant  coefficient  model  becomes  a  dynamic  coefficient  model 
as  the  model  coefficient  is  no  longer  taken  as  constant  but  allowed  to  vary  in 
space  and  time. 

Over  the  past  three  decades,  a  significant  amount  of  research  has  been 
devoted  to  LES  beginning  with  Lilly  [7]  who  adopted  LES  to  predict  flows  in 
the  field  of  Meteorology.  The  overwhelming  majority  of  this  research  has  been 
carried  out  with  spectral  or  structured  grid  finite  difference  methods.  However, 
recently  LES  has  been  extended  to  finite  element  methods  on  unstructured 
grids  (see  [4]).  This  extension  not  only  facilitates  simulation  of  flows  within  or 
around  complex  geometries,  but  also  allows  great  reductions  in  computational 
effort  through  the  ability  of  unstructured  grids  to  adapt  locally  to  resolve  fine- 
scale  flow  structures  in  one  flow  region  while  remaining  coarse  in  other  regions 
where  the  flow  structures  are  large. 

The  purpose  of  this  report  is  to  investigate  the  extension  of  LES  to  finite 
elements  by  focusing  on  the  intricacies  arising  from  such  a  task.  In  Section 
2,  we  present  classical  definitions  of  spatial  filters  along  with  their  relevant 
characteristics  for  LES.  In  Sections  3  and  4,  we  review  the  filtered  Navier-Stokes 
equations  (which  govern  the  large-scale  component  of  the  flow)  along  with  the 
dynamic  coefficient  Smagorinsky  model  required  to  represent  unknowns  in  the 
filtered  equations  generated  by  a  filtering  operation.  In  Section  5,  we  discuss 
discrete  filters  which  arise  when  the  modeled  filtered  Navier-Stokes  equations 
are  solved  using  finite  elements,  and  furthermore  we  introduce  new  filters  which 
make  use  of  the  basis  functions  underlying  the  numerical  method,  referred  to 
as  finite  dimensional  filters.  Finally,  in  Section  6  we  test  our  discrete  and  finite 
dimensional  filters  by  performing  large-eddy  simulations  of  decaying  isotropic 
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turbulence  and  comparing  the  results  with  experimental  data. 


2  Filtering 

In  the  large-eddy  simulation  (LES)  of  turbulent  flows,  the  larger  unsteady  tur¬ 
bulent  motions  (or  eddies)  are  directly  represented,  whereas  the  effects  of  the 
smaller-scale  motions  are  modeled.  There  are  four  conceptual  steps  in  LES: 

1.  A  filtering  operation  is  defined  to  decompose  the  velocity  u(x,  t)  into  the 
sum  of  a  filtered  (or  resolved)  component  ii(x,t)  and  a  residual  or  subgrid 
scale  (SGS)  component  u'{x,t).  The  filtered  component  represents  the 
motions  of  the  large  eddies. 

2.  The  evolution  equations  for  the  filtered  velocity  field  are  derived  from  the 
Navier-Stokes  equations.  These  equations  axe  of  the  same  form  as  the 
Navier-Stokes  equations,  with  the  momentum  equation  containing  an  un¬ 
known  residual  (SGS)  stress  tensor  that  arises  from  the  residual  motions. 

3.  Closure  is  obtained  by  modeling  the  residual  stress  tensor;  in  our  case, 
the  dynamic  coefficient  Smagorinsky  model  is  used. 

4.  The  modeled  filtered  equations  are  solved  numerically  for  u(x,t),  which 
provides  an  approximation  to  the  large-scale  motions  of  the  turbulent 
flow. 

The  general  filtering  operation,  introduced  to  LES  by  Leonard  [6],  is  defined  as 

u(x,t)  =  J  G{x,y)  u(y,t)  dy,  (1) 

where  integration  is  over  the  entire  flow  domain.  The  filter  kernel  G(x,y), 
located  at  y  =  x,  is  chosen  to  have  small  compact  support  in  y  thereby  reducing 
the  region  of  integration  to  be  much  smaller  than  the  flow  domain.  The  simplest 
filter  kernels  are  spatially  homogeneous,  symmetric  in  y  about  x  (leading  to 
the  expression  of  G  as  G(\x  -  y|)),  and  are  required  to  satisfy  the  following 
normalization  condition: 

J  G{r)  dr  =  1,  (2) 

guaranteeing  that  the  filtering  operation  preserves  constants.  Homogeneity 
simply  refers  to  the  fact  that  the  shape  of  the  kernel  remains  constant  as  we 
move  the  point  y  —  x. 

Symmetric  filter  kernels  can  be  defined  in  one  dimension.  The  extension 
to  the  three  dimension  vector  case  is  straight  forward.  Two  commonly  used 
homogeneous,  non-negative,  symmetric  filters  are  the  box  filter  and  the  Gaus¬ 
sian  filter.  Consider  a  random  function  f(x)  characterized  by  high  frequencies. 
With  the  box  filter,  the  filtered  function  j(x)  is  simply  the  average  of  f{x) 
in  the  interval  (x  -  A/2,x  +  A/2).  Thus  the  box  filter  kernel  takes  on  the 
value  of  1/A  over  the  previously  mentioned  interval  and  is  zero  elsewhere.  The 
Gaussian  filter  kernel  is  the  Gaussian  distribution  with  mean  zero  and  variance 
a2  =  A2/12.  This  value  of  the  variance  was  chosen  by  Leonard  so  as  to  match 
the  second  moment  of  the  Gaussian  filter  kernel  to  that  of  the  box  filter  kernel. 
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A  box  or  Gaussian  filtered  function  f(x)  follows  the  general  trends  of  f(x ), 
but  the  short  fluctuations  of  length  scales  smaller  than  A  have  been  damped. 
This  admits  a  natural  way  of  categorizing  filters  by  their  widths,  taken  as  A. 
Just  like  the  Gaussian  filter,  other  non-negative,  symmetric  filters  are  required 
to  match  their  second  moments  with  that  of  the  box  filter,  resulting  in  the 
following  general  definition  of  their  widths: 

/  roo  \  1/2 

A  =  ^12  J  r2G{r)  drj  .  (3) 

This  definition  is  extensively  discussed  by  Lund  [9].  The  previous  formula  will 
be  of  interest  to  us  when  defining  discrete  filtering  procedures  for  the  dynamic 
Smagorinsky  model. 


3  The  filtered  Navier-Stokes  equations 


The  incompressible  Navier-Stokes  equations  axe 


dm 


dxi 


duj 

dt 


dujUi 

dxj 


0 

d2Ui  1  dp 
dxjdxj  p dxi  ’ 


(4) 


where  the  first  equation  describes  conservation  of  mass  and  the  second  equation 
describes  conservation  of  momentum  in  the  i-th  direction.  As  usual,  Ui(x,t)  is 
the  component  of  the  velocity  field  in  the  i-th  direction,  p(x,t )  is  the  pressure, 
p  is  the  fluid  density  and  u  is  the  kinematic  viscosity.  Consider  the  appli¬ 
cation  of  any  homogeneous,  symmetric  filter  with  kernel  G  and  of  width  A. 
Consequently,  the  filtered  Navier-Stokes  equations  axe  rendered  as: 


du{ 


dxi 

dv,i  dujui 
dt  dxj 


0 

d2Ui  1  dp 
dxjdxj  p dxi 


(5) 


As  can  be  seen,  the  filtering  operation  commutes  with  differentiation,  which  is 
brought  about  by  the  homogeneity  of  the  kernels  in  question.  The  equations  in 
(5)  differ  from  the  Ne\  ier-Stokes  equations  because  the  filtered  product  uiu~j  is 
different  from  the  product  of  the  filtered  velocities  UiUj.  The  difference  is  the 
residual  or  subgrid  scale  stress  tensor  defined  by 


Tij  =  UiUj  —  UiUj. 


(6) 


Expressing  the  deviatoric  or  traceless  part  of  this  tensor  as 


rd  -  r  ■ 
Ti]  ~  t%3 


g  Tkk$ij  i 


(7) 


and  adding  the  trace  of  the  residual  stress  tensor  multiplied  by  the  density  to 
the  pressure  as 


P=P  +  P  gTfcfc, 


(8) 


5 


the  filtered  momentum  equation  in  (5)  can  be  rewritten  as 

dui  d'UjUi  _  y  d2Uj  _  drfj  _  1  dP 

dt  dxj  V  dijdxj  dxj  pdx{ 

If  rfj(x,t)  is  given  by  a  residual  stress  model,  the  filtered  continuity  equation 
(the  first  equation  in  (5))  and  the  filtered  momentum  equation  (9),  herein 
referred  to  as  the  filtered  Navier-Stokes  equations,  can  be  solved  to  determine 
u(x,  t)  and  P(x,t). 

In  the  case  of  compressible  flows,  the  compressible  Navier-Stokes  equations 
are  filtered  using  a  generalized  filtering  procedure  known  as  Favre  (or  density 
weighted)  filtering.  The  interested  reader  is  encouraged  to  consult  Moin  et  al 
[10]. 


4  The  dynamic  Smagorinsky  model 


The  filter  applied  to  the  Navier-Stokes  equations  is  meant  to  remove  the  small 
scales  of  motion.  Motion  is  converted  to  heat  via  viscous  dissipation  at  the 
smallest  of  these  small  scales.  Thus,  the  dissipation  of  motion  must  be  modeled 
due  to  the  absence  of  the  smallest  scales.  Fortunately,  by  Kolmogorov’s  first 
hypothesis  (see  [12]),  the  behavior  of  the  absent  smallest  (or  Kolmogorov)  scales 
is  universal,  and  as  a  result,  it  should  be  possible  to  construct  a  model  applicable 
to  all  types  of  flows.  A  simple  model  was  proposed  by  Smagorinsky  [15]  to 
account  for  the  dissipation  of  motion.  Smagorinsky  expressed  the  deviatoric 
portion  of  the  residual  stress  tensor  as 

rfj  =  -2(CA)2|S|%  (10) 


where  C  is  the  Smagorinsky  coefficient,  the  filtered  strain  rate  tensor  is  defined 
as 


a  l  ( duj  |  duj 

1J  2  ^  dx.j  dxi 


(11) 


and  the  norm  of  the  filtered  strain  rate  tensor  is  defined  as 


|5|  =  (2  SijSij)1'2.  (12) 

Note  that  (10)  is  consistent  in  the  sense  that  its  left  hand  side  and  its  right 
hand  are  traceless.  The  previous  model  is  valid  only  if  the  filter  width,  A, 
is  in  the  inertial  sub-range.  Scales  in  the  inertial  sub-range  are  larger  than 
the  Kolmogorov  scales,  yet  smaller  than  the  scales  which  contain  most  of  the 
energy.  Moreover,  the  inertial  sub-range  scales  are  universal  and,  unlike  the 
Kolmogorov  scales,  are  not  affected  by  molecular  viscosity.  The  reason  for  the 
restriction  on  A  is  that,  by  construction,  the  Smagorinsky  model  represents  all 
the  dissipation  of  motion  occurring  at  the  Kolmogorov  scales.  If  some  of  these 
scales  are  preserved  by  the  filter,  then  the  dissipation  of  motion  occurring  at 
these  scales  will  be  accounted  by  the  model  as  well.  This  is  clearly  incorrect. 
A  final  characteristic  of  the  Smagorinsky  model  is  that  it  is  not  suitable  to 
represent  the  possible  energy  transfers  from  the  removed  Kolmogorov  scales 


6 


to  the  larger  scales,  a  phenomenon  known  as  back  scatter .  The  reason  for 
such  deficiency  is  that  backscatter  occurs  on  time  scales  smaller  than  those 
represented  by  the  model.  Currently,  we  are  implementing  a  model  which  can 
exhibit  backscatter.  known  as  the  Bardina  mixed-scale  model  (see  [16]  and 

[19]))- 

One  of  the  problems  with  the  implementation  of  the  Smagorinsky  model 
is  that  the  appropriate  value  of  the  coefficient  C  is  different  in  different  flow 
regimes.  More  precisely,  it  is  zero  in  laminar  flow,  and  it  is  attenuated  near 
walls  compared  to  its  value  ( C  «  0.15)  in  high  Reynolds  number  free  turbulent 
flows.  To  alleviate  this  deficiency  of  the  constant  coefficient  model,  Germano 
[3]  developed  a  procedure  for  calculating  the  coefficient  locally.  Consider  the 
application  of  a  homogeneous,  symmetric,  secondary  or  test  filter  (with  kernel 
G  and  of  width  A  >  A  in  the  inertial  sub-range)  to  the  once  filtered  Navier- 
Stokes  equations.  As  expected,  the  continuity  and  momentum  equations  take 
the  form: 


dui 

dxi 

dui  duj7ui 
dt  +  dxi 


=  0 


=  v 


d2Uj 
dxj  d  x j 


1  dp 
pdx{ 


(13) 


The  residual  stresses  created  by  the  application  of  the  test  filter  are  expressed 


as 


Tij  —  UiU  j 


UiUj. 


(14) 


The  successive  application  of  two  homogeneous  filters  of  widths  A  and  A  re¬ 
spectively,  yields  a  homogeneous  filter  with  kernel  G  and  of  width  A  >  A 
(see  [12]  and  [17]).  (see  [16]  and  [18]).  Thus  the  Smagorinsky  model  for  the 
deviatoric  portion  of  the  residual  stresses  in  (14)  takes  shape  as 


rpd  _  rp 

~  1ij 


\ Tki kSij  =  -2(CA)2|5|4-, 


(15) 


where  §ij  and  \S\  are  defined  straightforwardly  based  on  u.  Moreover,  C  is  the 
same  constant  appearing  in  (10),  a  consequence  of  Kolmogorov’s  hypothesis 
mentioned  earlier. 

An  identity  due  to  Germano  is  obtained  by  applying  the  test  filter  to  (6) 
and  subtracting  the  result  from  (14): 


Lij  —  *7 ~ij  —  UiUj  UiU  j.  ( 1 6) 

The  significance  of  this  identity  is  that  (called  the  resolved  stress  tensor) 
is  known  in  terms  of  it,  whereas  and  are  not.  Taking  the  Smagorinsky 
coefficient  as  roughly  constant  in  the  neighborhood  where  the  filter  kernel  is 
non-zero,  the  relations  in  (10)  and  (15)  lead  to 


Lfj=Tfj-rfJ  =  2(CAfMij, 

(17) 

where  Lf-  is  the  deviatoric  portion  of  Lip 

^ij  —  Lij  —  —LkfcSij , 

(18) 

7 


and 


2 


Mij  =  |5|5j  •  -  x  |S|% 


(19) 


The  coefficient  (CA)2  cannot  be  chosen  to  match  the  six  independent  relations 
in  (17),  but  as  shown  by  Lilly  [8],  the  mean-square  error  is  minimized  by  the 
specification  that 


(CA)2 


1  _  1  MjjLjj 

2  MkiMki  2  Mfc/M*;/ 


(20) 


The  second  equality  in  the  previous  expression  arises  from  the  fact  that  M ij 
is  traceless.  Finally,  both  Lij  and  M{j  are  known  in  terms  of  u(x:t):  thereby 
completing  the  expression  in  (10). 

The  denominator  and  numerator  in  (20)  are  averaged  over  homogeneous 
directions  of  the  flow  (i.e.,.  directions  over  which  mean  quantities  of  the  flow  are 
constant)  resulting  in  a  time  dependent  Smagorinsky  coefficient  as  a  function 
of  the  directions  in  which  the  flow  is  not  homogeneous.  Averaging  is  an  ad 
hoc  operation  performed  as  a  way  to  avoid  instabilities.  In  some  instances,  the 
fluctuations  of  M{j  between  negative  and  positive  values  result  in  a  vanishing 
of  the  denominator  in  (20).  Furthermore,  a  negative  value  of  (C A)2  would 
give  rise  to  a  parasitic  transfer  of  energy  from  the  modeled  smallest  scales  to 
the  larger  scales.  We  refer  to  this  transfer  of  energy  as  parasitic  because,  as 
mentioned  earlier,  the  Smagorinsky  model  is  not  able  to  describe  such  transfers 
(back  scatter). 


5  Test  filtering 

Before  we  begin  our  discussion  on  test  filtering,  it  is  necessary  to  mention 
the  numerical  method  used  to  solve  the  filtered  Navier-Stokes  equations.  The 
spatial  discretization  in  our  large-eddy  simulations  is  brought  about  by  the 
Galerkin  approximation  to  the  weak  form  of  the  filtered  Navier-Stokes  equations 
augmented  by  the  Streamline  Upwind  Petrov- Galerkin  (SUPG)  stabilization, 
herein  referred  to  as  the  stabilized  Galerkin  approximation  (see  [4]  and  [18]). 
The  spatially  discretized  equations  are  advanced  in  time  using  the  generalized- a 
time  integrator  discussed  in  reference  [5].  The  stabilized  Galerkin  approxima¬ 
tion  involves  the  representation  of  the  flow  variables  as  a  linear  combination 
of  continuous  piecewise  polynomial  basis  functions  of  order  p.  In  turn,  these 
velocity  fields  are  used  to  calculate  the  dynamic  model,  which  requires  the  test 
filtering  of  numerous  flow  quantities  including  the  product  of  two  components 
of  the  velocity  (see  (16)).  For  example,  if  p  =  1  these  flow  quantities  are  prod¬ 
ucts  of  the  piecewise  tri-linear  linear  velocity  fields,  thus  they  will  be  piecewise 
tri-quadratics  or  higher. 


5.1  The  box  filter  as  the  test  filter 

Performing  the  filtering  integrations  using  the  box  filter  kernel  (of  width  A  —  2 h 
where  h  is  the  constant  mesh  spacing)  integrated  by  different  quadrature  rules 
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Nodal  point 


Quadrature  point 


Figure  1:  A  nodal  point  with  neighboring  quadrature  points  are  depicted.  Test  filtered 
quantities  are  calculated  at  the  nodes. 

admits  a  family  of  discrete  approximations  to  a  box  filtered  function  evalu¬ 
ated  at  a  particular  node  in  the  mesh.  In  one-dimension,  this  family  can  be 
represented  as 

fix o)  =  E  Wif(xi).  (21) 

i=-J 

Here  the  filtered  function  f{x)  is  evaluated  at  a  mesh  node  whose  spatial  lo¬ 
cation  is  denoted  as  xo  with  neighboring  quadrature  points  located  at  x\  = 
xo  +  I'i,  2:2  =  xo  +  L2,-  ••  ,xj  =  xo  +  Lj  to  the  right  and  x-i  =x0-Li,  x-2  = 
xo  ~  L2,  •  •  ■ ,  x-j  —  xq  -Lj  to  the  left,  where  the  {Li}  are  constants  determined 
by  the  quadrature  rule  and  the  mesh  spacing  h. 

The  {Li}  are  representable  in  the  form  Lj  =  Oj/i  where  0  <  a*  <  1.  Further¬ 
more,  the  weights  are  non-negative,  symmetric  (i.e.,  W{  =  W-i),  and  satisfy 
the  condition 

E  Wi  =  1.  (22) 

i=-J 

Therefore,  this  family  of  discrete  filters  preserves  constants.  An  example  of  a 
discrete  approximation  to  the  box  filter  is  given  in  Appendix  I.  The  family  of 
discrete  kernels  corresponding  to  the  family  of  discrete  filters  applied  in  (21) 
can  be  expressed  as 

J 

G(x  -  y)  =  h  E  Wi6(x  -y  +  a,-),  (23) 

i=-J 

where  =  ~oti  and  5(:r)  is  the  Dirac  delta  function.  The  previous  relation 
can  be  inserted  into  (3)  to  obtain  a  general  expression  for  the  filter  width: 

A  =  h  ^12  E  ■  (24) 

It  is  important  to  note  that  the  width  of  the  discrete  filters  is  proportional  to  the 
mesh  size  h.  Thus,  if  the  mesh  is  too  coarse,  the  width  of  these  discrete  filters 
might  be  too  wide  for  consistent  calculation  of  the  dynamic  model.  Recall  that 
the  width  of  the  test  filter  should  be  within  the  inertial  sub-range  discussed 
earlier. 


5.1.1  Discrete  filtering  in  multi-dimensions 

In  the  case  where  filtering  is  required  in  more  than  one  dimension  over  any 
topology,  a  box  filtered  function  evaluated  at  node  A  is  given  as 

f{xA)  - - 777-T  [  f{y)dy.  (25) 

meas{iiA)  JnA 

where  is  the  union  of  elements  which  share  node  A.  Here,  the  box  filter 
kernel  has  been  generalized  such  that  at  node  A  it  becomes 

G(xA,y)  =  |  meas(fi^)  lf  y  1S  1  A  (26) 

[  0  otherwise. 

In  the  case  where  filtering  is  performed  over  an  evenly  spaced  quadrilat¬ 
eral  or  hexahedral  rnesh,  quadrature  approximations  admit  the  sequential  ap¬ 
plication  of  the  one-dimensional  discrete  filters  in  each  of  the  three  principal 
directions.  For  example,  in  three  dimensions  we  have 

j  J  J 

f(xo,yo,zo)=  Z  £  £  WiWjWkfixuy^Zk),  (27) 

i=-J  j—  —  J  k——J 

where  the  {yj}  and  the  {zj}  are  spaced  in  the  same  manner  as  the  {xj}  de¬ 
fined  earlier.  Expanding  the  previous  sum,  the  filtered  function  at  the  node 
(zo,2/o,2o)  is  seen  to  be  given  by  a  weighted  combination  of  the  function  eval¬ 
uated  at  neighboring  quadrature  points.  The  weights  corresponding  to  each 
point  add  up  to  unity.  The  width  of  the  three  dimensional  filter  in  (27)  is 
defined  as  the  width  of  its  one  dimensional  counterpart. 

In  addition,  it  is  possible  to  sequentially  apply  filters  of  different  widths,  for 
which  the  width  of  the  resulting  filter  is  taken  often  as 

A  =  (A1A2A3)  ^  ,  (28) 

where  A*  is  the  width  of  the  filter  applied  in  the  i- th  direction.  In  the  case  of 
elements  with  high  aspect  ratios,  the  width  is  taken  as 

A  =  max  (Ax,  A2,  A3)  .  (29) 

For  a  detailed  discussion  on  when  it  is  appropriate  to  use  either  of  the  two 
previous  expressions  the  reader  is  suggested  to  consult  Deardoff  [2]  and  Scotti 
et  al  [13],  [14]. 

As  discussed  before,  the  test  filter  must  be  homogeneous  in  order  to  derive 
the  dynamic  model.  Specifically,  the  test  filter  must  be  able  to  commute  with 
differentiation.  Errors  induced  by  the  non-commutivity  between  differentiation 
and  discrete  filtering  operations  such  as  (21)  are  of  0(A2)  =  0(ft2),  which  is 
the  same  order  as  the  spatial  and  temporal  discretizations  in  the  flow  solver. 
A  detailed  discussion  on  the  commutation  errors  of  discrete  filters  is  provided 
by  Vasilyev  [16]. 

In  the  case  where  two  or  three  dimensional  filtering  is  performed  using  un¬ 
structured  (triangular,  tetrahedral,  non-uniform  quadrilateral,  or  non-uniform 
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hexahedral)  meshes,  the  generalized  box  kernel  may  not  be  symmetric;  conse¬ 
quently,  its  width  may  not  be  obtainable  from  the  formula  in  (3).  Furthermore, 
in  this  case  the  box  kernel  may  not  be  homogeneous,  thus  discrete  approxima¬ 
tions  may  induce  large  commutation  errors.  Later  in  this  report  we  will  show 
that  discrete  approximations  of  box  filter  on  uniform  tetrahedral  meshes  can 
work  well  in  large-eddy  simulations.  In  the  near  future  we  hope  to  show  that 
these  non-homogeneous  discrete  filtering  operations  introduce  negligible  com¬ 
mutation  errors. 

Furthermore,  in  the  case  of  unstructured  meshes,  the  quadrature  rules  used 
to  evaluate  the  filtering  operation  in  (25)  does  not  yield  filtered  functions  rep¬ 
resentable  in  the  form  of  (27).  However,  still  they  result  in  a  weighted  combi¬ 
nation  of  the  function  evaluated  at  the  points  surrounding  that  node  at  which 
the  filtered  function  is  evaluated.  By  construe!  ion,  the  weights  corresponding 
to  these  points  are  non-negative  and  satisfy  the  normalization  condition. 

5.2  Finite  dimensional  filters  as  test  filters 

We  hope  to  perform  large  eddy-simulations  using  higher  order  hierarchic  basis 
functions  in  addition  to  the  continuous  piecewise  linear  basis  functions  cur¬ 
rently  used.  The  purpose  for  using  higher  order  basis  functions  is  to  enable  the 
efficient,  accurate  representation  of  turbulent  flows  on  coarse  meshes.  If  we  use 
the  box  filter  as  the  test  filter  in  the  same  manner  it  was  used  in  the  previous 
sub-section,  the  resulting  filter  width  will  be  greater  than  the  inertial  sub-range 
scales  due  to  the  coarseness  of  the  meshes,  thereby  making  the  Smagorinsky 
model  invalid.  Thus,  we  need  to  define  a  special  class  of  filters  such  that  their 
widths  fall  on  the  inertial  sub-range  even  on  coarse  meshes.  For  this,  we  make 
use  of  ideas  discussed  by  Leonard  [6]  and  Pope  [11],  who  envision  a  filtered 
function  as  a  finite  dimensional  projection  of  an  infinite  dimensional  function. 

Once  again,  consider  a  random  function  /  (x)  characterized  by  high  frequen¬ 
cies.  The  motivation  for  using  finite  dimensional  projections  of  /( x)  as  filters 
is  that  they  can  represent  /  (x)  up  to  certain  scales.  More  precisely,  the  higher 
the  order  of  the  basis  functions  spanning  the  space  being  projected  onto,  the 
better  the  projected  function  can  represent  the  smaller  scales  of  f(x).  Such 
projections  can  be  expressed  as 


;=f>n<M-),  (30) 

n— 1 

where  the  members  of  the  set  {</>n}^=i  are  continuous  piecewise  polynomial 
basis  functions  and  the  set  of  coefficients  {on(x)}^=1  are  to  be  determined. 

5.2.1  Piecewise  polynomial  interpolations 

One  family  of  filters  can  be  obtained  by  selecting  the  set  of  coefficients  {a„  }^=1 
such  that  the  filtered  function  interpolates  the  original  function  through  N  pre¬ 
determined  points.  As  an  example  of  this  type  of  filters,  we  take  the  filtered 
function  f(x)  to  be  the  piecewise  linear  polynomial  interpolating  the  original 
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function  f(x)  through  the  mesh  nodes: 

f  =  *52  fn<t>n{x),  (31) 

n= 1 

where  fn  =  /(xn),  N  is  the  number  of  nodes  in  the  mesh,  and  the  set  {$n{%)}n=i 
consists  of  the  Lagrangian  continuous  piecewise  linear  basis  functions.  The  ker¬ 
nel  associated  to  this  filtering  operation  can  be  written  as 

N 

G(x ,  y)  =  ~  xn)<l>n{z)-  (32) 

71  =  1 

Notice  that  this  kernel  is  similar  to  those  in  (21),  except  that  instead  of  constant 
weights  it  has  variable  weights.  The  homogeneity  and  symmetry  of  the  kernel 
in  (32)  is  still  in  question.  However,  this  kernel  works  well  in  the  large-eddy 
simulation  of  decaying  isotropic  turbulence.  Similar  to  the  generalized  box  filter 
on  hexahedral  (or  quadrilateral)  elements,  we  hope  to  show  that  the  filtering 
operation  brought  about  by  (30)  introduces  negligible  commutation  errors. 

A  disadvantage  in  using  kernels  such  as  the  one  in  (32)  to  compute  the 
dynamic  model  coefficient  in  (20)  is  that  gradients  of  the  solution  at  mesh  nodes 
are  not  continuous  when  the  solution  is  piecewise  linear  or  any  other  piecewise 
polynomial.  However,  we  can  overcome  this  obstacle  by  assuming  that  the 
required  gradients  at  mesh  nodes  are  averages  over  surrounding  elements. 

5.2.2  Least-squares  (L2)  projections 

A  second  family  of  filters  can  be  obtained  by  determining  the  set  of  coefficients 
{«n}n"-i  hi  (30)  such  that  they  minimize  the  square  of  the  vT-norni  of  the 
difference  between  the  filtered  and  original  functions  over  the  domain  of  interest: 

R=I,Io  “  /(x))2  dx »  (33) 

where  L  is  the  distance  of  our  domain.  By  substituting  (30)  into  (33)  and 
differentiating  with  respect  to  the  {om},  it  is  seen  that  the  coefficients  satisfy 
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the  matrix  equation 


A' 


71=1 

(34) 

with 

1  fL 

Bjrin  =  Jq  0m(^)^n(^)  dx. 

(35) 

and 

1  fL 

Vm  =  T  dx. 

L  Jo 

(36) 

The  family  of  filters 

in  question  admits  the  following  family  of  kernels: 

N  N 

G{x,y)  =  Y,  E  Bmn<l>n{x)<l>m{y), 

71=1  771=1 

(37) 

where  B^n  denotes  the  m-n  entry  of  the  inverse  of  the  Gram  (mass)  matrix 
B.  It  is  important  to  note  that  the  minimization  of  the  integral  in  (33)  could 
be  performed  over  each  mesh  element,  instead  of  the  flow  domain.  In  this  case, 
the  resulting  filtered  field  would  be  discontinuous  across  elements. 

Although  the  widths  of  the  finite  dimensional  filters  proposed  are  unclear, 
it  should  reflect  the  polynomial  order  of  the  basis  functions  employed.  Clearly, 
higher  order  basis  functions  would  yield  smaller  filter  widths.  Thus,  unlike 
discrete  filter  widths,  which  are  solely  proportional  to  the  mesh  size,  finite 
dimensional  filter  widths  should  depend  on  the  chosen  basis  functions  as  well. 

6  LES  of  Decaying  Isotropic  Turbulence 

6.1  Decaying  Isotropic  Turbulence 

Before  moving  on  to  a  discussion  of  numerical  results,  let  us  take  some  time 
to  review  the  essential  aspects  of  our  test  problem.  Under  certain  conditions, 
large-scale  motions  can  become  turbulent.  In  other  words,  the  large-scale  mo¬ 
tions  become  unstable  and  break  into  smaller  scale  motions  which  take  energy 
from  the  larger  ones.  Energy  is  passed  down  to  such  small  scales  at  which  it  is 
dissipated  by  the  action  of  molecular  viscosity.  At  high  enough  Reynolds  num¬ 
bers,  the  small-scale  motions  cease  to  depend  on  the  nature  of  the  large-scale 
flow,  leading  to  the  universality  of  small-scale  motions  discussed  in  Section  3. 
Furthermore,  these  scales  lose  all  directional  orientation,  in  other  words,  they 
become  isotropic.  The  energy  contained  in  the  sub-inertial  scales  (discussed  in 
section  3)  is  characterized  by  what  is  usually  referred  to  as  the  five  -  thirds 
law.  In  other  words,  the  energy  at  these  scales  behaves  as  kT  ,  where  kr 
represents  the  inverse  of  the  size  of  the  scales. 

In  the  up-coming  sub-sections  we  try  to  simulate  a  flow  which  is  nearly 
isotropic  at  all  scales.  Our  results  are  compared  to  the  experimental  data 
of  Comte-Bellot  and  Corrsin  [1],  who  tried  to  represent  an  infinite  space  of 
isotropic  motions  decaying  in  time  because  of  a  lack  of  kinetic  energy  produc¬ 
tion  (in  the  absence  of  shear  flows)  to  balance  the  viscous  dissipation.  They 
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accomplished  this  by  obtaining  a  turbulence  field  behind  a  regular  grid  span¬ 
ning  a  steady,  uniform  duct  flow.  By  moving  at  the  speed  of  the  mean  flow 
behind  the  grid,  they  correctly  surmised  that  an  observer  would  see  something 
like  true  isotropic  turbulence  evolving  in  time. 


6.2  Numerical  observations 

As  can  be  seen  in  the  dynamic  Smagorinsky  model,  the  filtered  velocity  field 
u  depends  on  the  ratio  A/ A.  The  primary  filter  kernel,  G ,  appears  implicitly 
in  the  equations,  while  the  kernel  G  appears  explicitly  when  computing  the 
dynamic  model.  Essentially,  the  primary  filter  kernel  as  well  as  the  kernel  G 
are  arbitrary,  yet  their  widths  appear  in  the  Smagorinsky  model.  Customarily, 
A  is  given  by  the  mesh  spacing  in  a  hexahedral  mesh  simulation: 

A  =  (fci/i2/i3)1/3  (38) 

or,  in  the  case  of  elements  with  high  aspect  ratios, 

A  =  max  (h\,  /i2,  ^3)  5  (39) 

where  hi  is  the  mesh  spacing  in  the  i-th  direction.  The  idea  behind  these 
choices  is  that  they  characterize  the  scales  at  the  resolution  threshold  of  the 
mesh.  In  reality,  when  the  filtered  equations  are  discretized,  A  as  well  as 
A  become  dependent  upon  the  numerical  method  of  choice  coupled  with  the 
effectiveness  (or  ineffectiveness)  of  the  dynamic  model,  the  mesh  spacing,  and 
the  polynomial  order  of  the  basis  functions  used  in  the  stabilized  Galerkin 
approximation,  which  together  act  as  the  primary  filter.  In  other  words,  A  in 
(38)  is  perturbed  in  the  form 

A  =  e(/ii/i2/i3)1/3,  (40) 


where  e  ~  1  is  a  constant  depending  on  the  numerical  method,  the  model,  the 
mesh,  and  the  polynomial  order  of  the  basis  functions  underlying  the  method. 
Thus,  we  are  left  with  a  difficult  question:  How  can  we  evaluate  A  and  A 
correctly?  More  precisely,  how  can  we  compute  the  ratio  A/ A  correctly?  In 
section  4,  we  showed  how  to  calculate  A  for  discrete  test  filters  on  a  uniform 
hexahedral  mesh.  Here  we  make  use  of  this  information  in  assuming  that  for  a 
hexahedral  mesh 

A_  /A1a2A3\1/3 
P  A  K{  h1h2h3  j  ’ 


(41) 


where  k  is  a  positive  constant  and  /3  is  the  filter  width  ratio  (FWR).  The 
parameter  k  serves  to  characterize  the  unknown  filter  width  A. 

In  the  decay  of  isotropic  turbulence,  the  scales  of  motion  are  isotropic  and 
vary  over  the  same  range  everywhere  in  the  domain.  Periodicity  can  be  as¬ 
sumed  in  the  three  principal  directions  as  long  as  the  two-point  spatial  auto¬ 
correlations  of  velocities  over  distances  equivalent  to  the  domain  length  vanish. 
Given  these  characteristics  of  the  flow  in  question,  we  can  ascertain  that  for 


14 


an  evenly  spaced  hexahedral  mesh  with  no  p-refinement  the  numerical  pro¬ 
cedure  resolves  turbulent  motions  larger  than  a  constant  scale  A.  (Otherwise, 
if  we  were  to  perform  p-refinement  or  h-refinement  in  some  regions  of  this 
flow,  the  numerical  procedure  would  resolve  up  to  scales  smaller  than  in  other 
regions,  yielding  a  non-constant  A.)  Moreover,  A  is  constant  as  well  because 
we  are  using  an  evenly  spaced  mesh.  Thus,  A,  which  depends  on  A  and  A,  is 
taken  as  constant.  Now  it  should  be  clear  that  for  decaying  isotropic  turbulence 
simulated  on  an  evenly  spaced  hexahedral  mesh  with  no  p-refinement ,  the  ex¬ 
pression  in  (41)  is  consistent  in  the  sense  that  both  of  its  sides  are  constants 
throughout  the  flow  domain. 

To  verify  the  validity  of  (41),  we  performed  the  simulation  with  two  different 
test  filters  on  a  fixed  evenly  spaced  hexahedral  mesh  using  continuous  piecewise 
linear  basis  functions.  For  each  case,  the  filter  width  ratio  was  expressed  as: 


Ai  Ax 
=  «- 


A2  A2 

and  — . 

A  h 


(42) 


A  "  h 

Here  the  subscript  on  the  filter  width  denotes  the  first  or  second  simula¬ 
tion.  Note  that  no  distinction  is  made  between  the  filters  applied  in  dif¬ 
ferent  directions  because  we  are  dealing  with  an  evenly  spaced  mesh  (i.e.,. 
h  =  h\  =  /i2  —  ^3)-  Furthermore,  no  distinction  is  made  for  A  between  the 
first  and  second  simulations  arguing  that  because  the  uniform  mesh  is  fixed,  A 
remains  constant  as  was  discussed  earlier.  However,  one  might  also  argue  that 
because  the  test  filter  changes  the  numerical  method  changes  as  well,  conse¬ 
quently,  A  changes.  Later  we  will  see  that  this  is  not  the  case  as  long  as  the  test 
filter  widths,  Ai  and  A2,  are  consistently  calculated  using  (24).  Thus,  because 
the  constant  k  characterizes  A,  it  remains  constant  as  well  from  simulation  to 
simulation.  In  summary,  we  will  see  that  both  A  and  k  are  independent  of  the 
test  filter  used  as  long  as  the  test  filter  width  is  accurately  computed. 

When  the  simulation  is  performed  with  a  uniform  tetrahedral  mesh,  the 
assumption  made  in  (41)  is  not  feasible  because  A  is  not  well-defined.  However, 
we  can  measure  the  constant  A 

fi  =  T  <43> 


for  which  the  simulation  perforins  well  with  respect  to  experimental  data.  Note 
that  the  best  value  of  /3  is  particular  to  the  test  filter.  Furthermore,  in  setting 
(43)  as  a  constant  we  have  assumed  that  for  an  uniform  tetrahedral  mesh  with 
no  p-refinement,  A  is  constant  throughout  the  flow  domain,  as  was  discussed 
earlier. 

In  the  case  of  more  complex  flows,  where  the  scales  of  motion  might  vary 
from  region  to  region,  h-refinement  and  p-ref inement  are  often  employed  to 
improve  local  resolution.  If  these  refinements  are  made,  no  longer  are  A,  A, 
and  A  constants  over  the  domain;  instead,  they  reflect  the  local  mesh  spacing 
and  local  polynomial  order.  However,  if  we  were  to  refine  such  that  both  sides 
of  (41)  or  the  left  hand  hand  side  of  (43)  remain  constants,  the  consistency 
of  these  expressions  should  hold,  at  least  approximately.  Their  validity  should 
be  verified.  Furthermore,  A  and  A  should  be  varied  slowly  enough  such  that 
they  can  be  approximated  as  locally  constant  in  the  derivation  of  the  dynamic 


model. 
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Figure  3:  LES  of  decaying  isotropic  turbulence  performed  with  hexahedral  elements. 
FWR  (filter  width  ratio)  =  /3. 

6.3  Numerical  results 

Two  sets  of  simulations  of  decaying  isotropic  turbulence  were  performed  using 
an  evenly  spaced  linear  hexahedral  mesh  discretizing  a  cubic  domain  with  sides 
of  length  27t  and  subjected  to  periodic  boundary  conditions  in  the  three  princi¬ 
pal  directions.  The  two  sets  are  distinguished  solely  by  a  change  in  the  discrete 
test  filter.  In  the  first  set,  the  test  filter  is  taken  as  the  box  filter  approximated 
by  one-point  quadrature.  In  the  second  set,  the  test  filter  is  taken  as  the  box 
filter  approximated  by  eight-point  quadrature. 

The  length  of  the  mesh  spacings  was  taken  as  h  =  2n/32,  where  32  is  the 
number  of  space  intervals  in  each  direction.  Furthermore,  as  mentioned  earlier, 
the  results  are  compared  to  experimental  data.  However,  since  the  results 
are  in  terms  of  primary  filtered  variables  ('.e..  u),  the  experimental  data  used 
for  comparison  should  also  be  filtered.  TL«  filtering  of  the  experimental  data 
should  be  of  width  A,  however,  since  this  value  is  not  available  to  us,  we  picked 
it  to  be  0(h),  specifically,  A  =  h. 

6.3.1  Box  filter  on  hexes  approximated  by  one-point  quadra¬ 
ture 

The  first  set  of  simulations  was  performed  with  the  test,  box  filter  approximated 
by  a  one-point  quadrature  rule  yielding  a  discrete  filter  of  width  A  =  \/3  h. 

First,  k  in  (41)  was  set  to  unity,  yielding  a  filter  width  ratio  (3  =  \/3-  Then 
k  was  set  to  any  number,  here  taken  as  ^/9/12,  yielding  a  filter  width  ratio 
/3  =  ,/9/4  =  \/2.25.  The  results  of  this  set  of  simulations  are  plotted  in  Figure 
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Energy  spectras  at  t=98  using  hexes  and  8-pt.  quad,  on  the  filter  and  the  flow 


Figure  4:  LES  of  decaying  isotropic  turbulence  performed  with  hexahedral  elements. 
3. 

6.3.2  Box  filter  on  hexes  approximated  by  eight-point  quadra¬ 
ture 

The  second  set  of  simulations  was  performed  using  the  test  (or  box)  filter  ap¬ 
proximated  by  an  eight-point  quadrature  rule,  yielding  a  discrete  filter  of  width 
A  =  2 h.  First,  k  was  set  to  unity,  yielding  a  filter  width  ratio  3  —  \/4-_Then, 
once  again,  k  was  reset  to  x/9/12,  yielding  a  filter  width  ratio  3  =  \/3-  The 
results  of  this  set  of  simulations  are  plotted  in  Figure  4. 

Looking  at  Figures  3  and  4,  we  see  that  results  do  not  vary  from  simulation 
to  simulation  as  long  as  the  test  filter  widths  axe  consistently  calculated  using 
(24).  We  ascertain  that  the  test  filter  does  not  have  an  impact  on  the  results 
as  long  as  its  width  is  accurately  represented.  Furthermore,  we  can  say  that 
A  remained  constant  from  simulation  to  simulation  as  was  assumed  before  and 
that  the  parameter  k  strictly  characterizes  the  unknown  primary  filter  width  A 
independent  of  the  test  filter  used.  In  addition,  the  ratio  3  scales  as  0(A/A) 
for  this  problem  and  the  parameter  k  can  be  adjusted  so  as  to  obtain  preferable 
results. 

6.3.3  LES  on  tetrahedral  elements 

The  simulation  of  decaying  isotropic  turbulence  was  also  performed  with  uni¬ 
form  tetrahedral  elements  discretizing  the  domain  discussed  in  the  previous 
paragraph.  As  mentioned  earlier,  a  proper  value  of  3  in  this  case  is  not  within 
our  means,  thus  we  chose  to  keep  this  value  fixed  at  \/3  for  all  simulations.  How- 
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Energy  spectras  at  t=98  using  tets  with  FWR=sqrt(3) 


Figure  5:  LES  of  decaying  isotropic  turbulence  performed  with  tetrahedral  elements 
and  ft  =  \/3. 


ever,  we  varied  the  quadrature  rule  used  to  approximate  the  box  filter  and  the 
quadrature  rule  used  to  integrate  the  discretized  flow  equations  (which  alters 
the  numerical  method).  Although  not  shown  in  the  figures,  a  slight  variation 
of  ft  affects  the  energy  spectra  shown  in  Figures  5  and  6,  thus  a  value  of  (3  can 
be  obtained  such  that  the  simulation  nearly  matches  experimental  results.  In 
the  future  we  hope  to  arrive  at  a  concrete  way  of  calculating  the  best  choice 
for  p.  Of  greater  importance  is  that  despite  the  possible  non-homogeneity  and 
non-symmetry  of  the  generalized  box  filter  on  tetrahedral  elements,  discrete 
approximations  of  said  filter  are  shown  to  work  just  as  well  as  discrete  ap¬ 
proximations  of  the  generalized  box  filter  on  hexahedral  elements,  for  which 
homogeneity  and  symmetry  are  well  established. 

6.3.4  An  approximate  least-squares  projection 

In  Figure  7  we  show  results  of  the  simulation  performed  with  the  least-squares 
projection  filter  whose  kernel  appears  in  (37).  This  filter  operator  was  ap¬ 
proximated  by  lumping  of  the  Gram  (or  mass  matrix)  and  using  one-point 
quadrature.  Here  we  have  chosen  (3  =  \/3,  for  which  the  simulation  performed 
well  with  hexahedral  elements.  Choosing  (3  =  \/4  produces  even  better  results. 
However,  in  the  future  we  would  like  to  establish  a  rigorous  way  of  how  to 
choose  the  appropriate  value  of  (3.  Although  not  shown,  the  simulation  also 
performed  well  with  tetrahedral  elements. 

Note  that  in  this  simulation  filtered  functions  where  evaluated  in  the  interior 
of  the  elements,  unlike  the  simulations  performed  with  the  discrete  box  filter  in 
which  the  filtered  functions  were  evaluated  at  the  mesh  nodes.  The  reason  for 
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Figure  6:  LES  of  decaying  isotropic  turbulence  performed  with  tetrahedral  elements 
and  P  —  y/3. 


this  variation  is  that  the  least-squares  filter  operator  under  the  lumping  and 
one-point  quadrature  approximations  reduces  to  the  box  kernel  approximated 
by  one-point  quadrature  when  both  axe  used  at  the  mesh  nodes. 

6-4  Initial  conditions  and  implementational  details 

Experimental  data  is  available  from  [1],  however,  it  is  not  enough  to  create  a 
fully  turbulent  initial  condition  for  the  simulations.  From  the  data  we  were 
able  to  obtain  the  correct  amplitudes  of  the  flow  velocities  in  wavenumber 
space  at  two  dimensionless  time  stations,  namely  t  =  42  and  t  =  98.  As  an 
approximate  initial  condition  we  used  the  experimental  velocity  amplitudes  at 
t  =  42,  at  discrete  points  in  wavenumber  space  corresponding  to  discrete  points 
in  real  space;  the  corresponding  phases  where  assigned  random  numbers.  These 
velocities  where  then  used  to  compute  a  pressure  field  satisfying  the  Poisson 
equation  for  pressure  in  wavenumber  space.  The  inverse  Fourier  transform 
was  used  to  obtain  real  initial  conditions.  However,  these  initial  conditions 
were  not  appropriate  due  to  the  random  phases  assigned  to  the  velocities.  The 
simulation  was  run  for  a  couple  of  hundreds  time  steps.  The  resulting  velocities 
where  then  transformed  back  to  wavenumber  space  and  their  phases  extracted 
and  re-assigned  as  the  phases  of  the  experimental  velocity  amplitudes  at  t  =  42. 
The  simulation  was  re-run  until  iterations  of  the  procedure  previously  outlined 
produced  the  same  results,  signaling  that  an  appropriate  (fully  turbulent)  initial 
condition  had  been  obtained. 

As  was  mentioned  earlier,  the  finite  element  solver  employed  uses  the  Stream- 
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Energy  spectra  at  t=98  using  hexes  and  the  approx,  l^-projection  filter 


Experimental  Data 
LES  w /  FWR=sqrt(3) 


Figure  7:  LES  of  decaying  isotropic  turbulence  using  an  approximation  of  the  least- 
squares  (L2)  projection  as  the  test  filter  and  f3  =  va¬ 


line  Upwind  Petrov-Galerkin  (SUPG)  method.  The  large  eddy-simulations  in 
this  report  were  performed  using  incompressible  (see  [18])  and  compressible 
(see  [4])  versions  of  the  solver.  The  compressible  simulations  were  performed 
with  a  Mach  number  much  less  than  unity  in  order  for  the  flow  to  approach 
incompressibility.  Furthermore,  the  compressible  solver  was  used  in  its  incom¬ 
pressible  limit.  Results  from  both  of  the  cases  using  the  compressible  solver 
matched  the  results  of  the  incompressible  solver. 

For  clarity  issues,  in  this  report  the  dynamic  Smagorinsky  model  was  derived 
in  its  incompressible  form.  For  compressible  flows  the  Smagorinsky  model  varies 
slightly.  The  interested  reader  is  encouraged  to  consult  [10]  for  a  derivation  of 
the  compressible  version  of  the  model.  This  compressible  version  reduces  to 
the  incompressible  form  by  setting  density  equal  to  a  constant. 


6.5  Conclusions 

Large-eddy  simulation  using  a  dynamic  residual  stress  model  proved  to  be  suc¬ 
cessful  when  extended  to  unstructured  meshes  using  a  generalized  box  filter  as 
the  test  filter.  Simulation  results  were  shown  to  be  independent  of  the  test  filter 
provided  the  test  filter  width  is  consistently  represented.  We  hope  to  further 
validate  the  role  of  the  parameter  k,  characterizing  the  unknown  primary  filter 
width,  by  performing  the  simulation  on  finer  meshes.  Furthermore,  possible 
alternative  filters  to  the  generalized  box  filter  were  introduced.  These  finite 
dimensional  filters,  based  on  continuous  piecewise  polynomial  approximations 
and  least-squares  projections,  would  be  preferred  for  higher  order  large-eddy 
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Figure  8:  Schematic  for  a  discrete  approximation  to  a  box  filtered  function. 

simulations  on  coarse  meshes.  In  the  near  future  we  hope  to  estimate  the 
homogeneity  of  the  generalized  box  filter  and  the  finite  dimensional  filters. 

7  Appendix  I 

Here  we  will  show  some  calculations  done  when  computing  a  box  filtered  func¬ 
tion  using  the  well  known  2-point  Gaussian  quadrature  rule.  Consider  the 
schematic  in  figure  8  as  we  set  out  to  obtain  the  filtered  function  /  evaluated 
at  node  A  or  0.  Here  xa  ~  %A-\  =  xa+ i  —  %a  =  h,  L\  =  h  ■  a\  =  h  ■  0.211  and 
1,2  =  h  •  d2  =  h  ■  0.788  Thus,  we  have 

2 

f(xA)  =  /(so)  =  E  Wif{xi),  (44) 

i— — 2 

where  W*  =  1/4  for  all  i.  This  discrete  filter  yields  a  filter  kernel  of  the  form 

2 

G{x-y)  =  h  E  Wi6{x-v  +  ai).  (45) 

i=— 2 

Inserting  the  values  for  Wj  and  a ;  with  J  =  2  into  the  expression  in  (24)  we 
find  that  the  width  of  the  filter  kernel  in  (45)  is  A  =  2 h. 
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